Reentrant phase transition in a predator prey model 
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We numerically investigate the six-species predator-prey game in complex networks as well as in 
d- dimensional regular hypercubic lattices with d = 1,2, •■■ ,6. The food-web topology of the six 
species contains two directed loops, each of which is composed of cyclically predating three species. 
As the mutation rate is lowered below the well-defined phase transition point, the Z2 symmetry 
related with the interchange of the two loops is spontaneously broken, and it has been known that 
the system develops the defensive alliance in which three cyclically predating species defend each 
other against the invasion of other species. In the Watts-Strogatz small-world network structure 
characterized by the rewiring probability a, the phase diagram shows the reentrant behavior as a is 
varied, indicating a twofold role of the shortcuts. In d-dimensional regular hypercubic lattices, the 
system also exhibits the reentrant phase transition as d is increased. We identify universality class 
of the phase transition and discuss the proper mean-field limit of the system. 



PACS numbers: 89.75.Hc, 89.75.Fb, 68.35.Rh, 87.23. Cc 



I. INTRODUCTION 



Recently, noncquilibrium dynamical systems as well as 
standard equilibrium statistical mechanical model sys- 
tems have been intensively studied in the interaction 
structure of complex networks [l|. Effects of shortcuts 
in the Watts-Strogatz (WS) networks [3| on collective 
dynamic behaviors have been drawn much attention, re- 
vealing the close interplay between the structural and the 
dynamical properties. For networks of highly heteroge- 
neous degree distributions such as the Barabasi-Albert 
(BA) scale- free network Q , various aspects of collective 
behaviors have also been studied. Lots of existing stud- 
ies have been performed in the framework of the phase 
transition and critical behavior [l| , and complex network 
topology in many cases gives rise to the mean-field uni- 
versality class. It has been agreed that the shortcuts 
and the hub vertices play important roles, increasing 
the effective dimensionality of the WS and the BA net- 
works, respectively, yielding the critical behavior beyond 
upper critical dimension. In contrast, although various 
noncquilibrium models like game theoretic models, epi- 
demic spread models, the voter model, and the contact 
process, have been actively studied in the complex net- 
work research area [l| , generic understanding of how the 
underlying network topology affects the nonequilibrium 
phase transition is still lacking. 

In population genetics, the so-called Lotka-Volterra 
model has often been studied. In the point of view of 
statistical physics, the neglect of the spatial density fluc- 
tuation in the original Lotka-Volterra model corresponds 
to a mean-field approximation, which cannot be justi- 
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fied in general, since every living organism inevitably 
lives in a finite-dimensional space with a limited range 
of interactions. Accordingly, population genetics mod- 
els allowing spatial inhomogeneity in species densities 
arc desirable, which can be simply realized by putting 
predators and preys on regular lattices in two- and three- 
dimensions (2D and 3D). Indeed, a cyclically interacting 
three-species predator-prey model has been studied both 
experimentally and numerically in Ref. y to uncover how 
the biodiversity of the system can be maintained. Sim- 
ilar three-species predator-prey model has also been in- 
vestigated not in a regular lattice structure but in the 
small- world interaction structure [a|- More coniplicated 
predator-prey model has been suggested later [a, 0, @1, 
with six or nine species interacting with each other in a 
given food-web structure. In these studies it has been 
found that a subgroup of cyclically predating species is 
spontaneously formed and species within the subgroup 
protect the member species from the attacks by other 
species outside of the subgroup. The spontaneous for- 
mation of such a defensive alliance is well captured by 
the statistical mechanical approach and the research fo- 
cus has been put on the existence and the nature of the 
phase transition of the spontaneous formation of alliance 
as the mutation rate is varied. 

In the present work, we study the six-species predator- 
prey game on various spatial interaction structures, in the 
context of the phase transition and the critical behavior 
of a noncquilibrium statistical physics model in complex 
networks. The food web under consideration (shown in 
Fig. [T]) has first been introduced in Ref. l2 and later stud- 
ied in Ref. [1. Specifically, Szabo et al. [3 have shown 
that the predator-prey game played on a 2D regular lat- 
tice exhibits a phase transition of the 2D Ising universal- 
ity as the mutation rate is changed, and also identified 
the relevant order parameter detecting the Z2 symmetry 
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FIG. 1: (Color online) Food web of the six-species predator- 
prey model. Each species has two predators, two preys, and 
two noninteracting neutral species. The two groups, called 
defensive alliances (alliance I and II), can be formed. 



breaking of the defensive alliance of three cyclically pre- 
dating species (see Fig.[T]). In what follows, we will refer 
to this model as the defensive alliance process (DAP). 
In Ref. 8 the DAP has been studied on the WS network 
structure and the effects played by the two different types 
of randomness, i.e., the temporal randomness induced by 
the mutation, and the structural randomness introduced 
by the shortcuts in the WS network, have been investi- 
gated. It has also been observed that the DAP in the 
WS network shows a discontinuous phase transition at 
nonzero rewiring probability different from the continu- 
ous transition in 2D regular lattice Q. We extend in 
this work the study in Ref. [a for the WS networks and 
construct the complete phase diagram in the plane of the 
two different randomness. To investigate the role of the 
structural inhomogeneity, the nature of the phase tran- 
sition in the regular d-D hypercubic lattices is also nu- 
merically studied. Since the mean-field theory developed 
in Ref. |8| does not predict a nontrivial critical point, we 
perform the cluster mean-field (CMP) approximation for 
2D model with the hope for the better understanding as 
to the nature of the transition in higher dimensions. 

The present paper is organized as follows: Section HIl 
presents our results of the phase diagram for the DAP 
on the WS networks as well as on the d-dimensioual hy- 
percubic lattices with d — 1, 2, ■ • • ,6 and discusses the 
reentrancc transition and the nature of the phase transi- 
tion. In Sec. mil the cluster mean-field theory with the 
cluster size 2 is developed for the 2D defensive alliance 
model. Finally, we summarize our results in Sec. IIVI 



II. DEFENSIVE ALLIANCE PROCESS 

We start from the description of how the DAP is im- 
plemented on general networks. The algorithm for the 
simulations is as follows: Initially at time t = 0, six 
species are equally distributed on vertices of the given 
network structure, with the density Cs of the species s 



given by Cs{t = 0) = 1/6 for all species (s = 0, 1, • • ■ ,5). 
At each time step, one vertex is chosen at random and 
with the probability P the species at the vertex is mu- 
tated to one of its predators. Otherwise, with probability 
1 — P, the species at the vertex plays the predator-prey 
game according to the rules depicted in Fig. [1] with one of 
its randomly selected neighbors, and the winner between 
the two occupies the loser's vertex. If the two species 
are neutral, i.e., if no arrow connects the two in Fig. [1] 
the game will end in a draw and nothing happens. The 
above procedure is repeated until the system approaches 
the steady state. 

Only for convenience, we define the parameter ^ as [8| 



M = ln(l/P), 



(1) 



which is a decreasing function of P. As P becomes larger, 
the temporal randomness becomes stronger. In other 
words, fi resembles an inverse temperature of equilibrium 
systems and 1/ ^ can be interpreted as an effective tem- 
perature. 

In order to detect the alliance breaking transition we 
measure the order parameter m (we also call it the mag- 
netization in analogy to the ferromagnetic Ising model) 
defined by ^M 



(|(C0 + C2+ Ci) - (Ci -f C3 -I- Cs)!), 



(2) 



where (...) is the time average after the steady state is 
achieved. In this work, a sufficiently long equilibration 
time (20 000 Monte-Carlo steps) is taken. As P becomes 
larger toward unity, all species are equally probable and 
m{P — > 1) = 0, while as P approaches zero, the sponta- 
neous development of defensive alliances gives us m « 1 
both for the alliances I and II (see Fig. [1]), indicating 
the possibility of a phase transition at nontrivial critical 
point lie- 

In simulating DAP, we had a numerical difficulty espe- 
cially when the mutation probability P is very small. In 
such cases, the system often spends very long time before 
achieving the steady state, which we try to avoid through 
the use of the simulated annealing technique in statistical 
physics by lowering P slowly starting from a high value of 
P. It is also observed that if the system size is not suffi- 
ciently large, the population becomes monomorphic and 
all vertices are occupied by a single species before mu- 
tation acts. The rare mutation can produce a predator 
of the species, but the small system at the low mutation 
probability will again be monomorphic unless a predator 
of the new born predator is generated by another muta- 
tion. The results presented below are for sufficiently large 
systems where the polymorphic population is attained. 



A. DAP in the WS networks 

We first construct WS networks for the DAP as fol- 
lows [2|: (i) ID and 2D regular lattices are first built. 
For ID each vertex has connections to its four neigh- 
bors (nearest and next-nearest neighbors), whereas for 
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FIG. 2: (Color online) The order parameter m detecting the 
alliance breaking transition versus the mutation parameter 
fj, = ln(l/P) with the mutation probability P is shown for 
the network WSi, constructed from the ID regular lattice via 
rewiring, at the rewiring probability (a) a = 0.02 and (b) 
a = 0.3. As /i is increased, i.e., as the mutation probability P 
is decreased, a spontaneous formation of the defensive alliance 
occurs. The transition point ^c is roughly estimated from the 
comparison of different sizes (iV = 60 000, 120 000, 240 000): 
fic = 7.4 ± 0.1 for (a) a = 0.02 and /ic = 6.8 ± 0.1 for (b) 
a = 0.3. 



2D we assume only four nearest neighbor connections. 
(ii) With the rewiring probability a one end of each local 
Unk is moved to a randomly chosen other lattice point. 
Throughout this paper, we call the resulting networks as 
WSi and WS2, respectively, in order to indicate the di- 
mensionality 1 and 2 of initial regular lattices from which 
the WS network is built as described above. However, 
the subscripts 1 and 2 in WSi and WS2 should not be 
interpreted as the dimensionality of the resulting small- 
world networks. The structure of the network is varied 
with a from a regular network (a = 0) to a fully random 
network {a ~ 1). 

Except for the case of a = 0, corresponding to locally 
connected ID and 2D regular lattices, the phase transi- 
tion is found to be of the discontinuous nature as was 
already found in Ref S When /x > /Zc, the order param- 
eter m saturates toward a finite value from below as the 
system size N is increased, while at n ^ fXc, it decreases 
with N, which implies that in the thermodynamic limit 
of A^ -^ 00, the order parameter changes abruptly, sig- 
nalling a discontinuous phase transition (sec Ref. (SJ for 
more detailed discussion). For examples, m versus fi for 
the WSi network is shown in Fig. [2] for a = (a) 0.02 and 
(b) 0.3, with the estimation fi^ = 7.4(1) and Hc = 6.8(1), 
respectively. We use this finite-size behavior to locate ^c 
as above, which is then used to construct the full phase 
diagrams in Fig. [3l 

It is very interesting that the phase diagrams in Fig. [3] 
exhibit reentrant behaviors as a is increased. This indi- 
cates that the role of the random shortcuts is somehow 
twofold, facilitating the defensive alliance at small a and 
then suppressing the formation of the alliance at large a. 
In detail. Fig. [3jc) displays m versus the rewiring prob- 
abihty a at fixed fi ^ 6.0 (l//i w 0.17) for the WS2 net- 
work, following the vertical arrow indicated in Fig.[3Jb). 
It is shown clearly that as a increases the system starts 



Disordered. 



(a) WSi 



Ordered 



Disordered. 
Jb)WS2 - 



Ordered 



0.1 0.2 0.3 

1/H 



0.1 0.2 0.3 

1/H 




FIG. 3: (Color online) Phase diagrams for the DAP in the 
plane of the mutation parameter l//i[= l/ln(l/P)] with the 
mutation rate P and the rewiring probability a for (a) WSi 
and (b) WS2 networks. As the rewiring probability a is 
changed, both (a) and (b) show reentrant phase transitions. 
A discontinuous phase transition is observed at any nonzero 
value of a. At a = 0, WSi and WS2 correspond to the regu- 
lar ID and 2D lattices, (c) m versus a for the WS2 networks 
of the sizes TV = 256 x 256 and 512 x 512 with /i set to 6.0, 
following the vertical arrow in (b). The reentrant behavior is 
again seen very clearly. 



from a disordered phase with a very small m and en- 
ters an ordered phase and then finally leaves back to a 
disordered phase. We also check the finite-size effect by 
comparing m for two different sizes N = 256 x 256 and 
512 X 512 in Fig.[31[c); the existence of ordered phase in 
the intermediate region of a is shown to be not a finite- 
size artifact. As a is increased from zero, more short- 
cuts make the system more strongly correlated in the 
sense that the change of the dynamic state of one vertex 
can affect more vertices due to the small- world effect [2[ . 
Consequently, we believe that the first increase of m for 
small a can be attributed to the strengthened correlation 
due to more shortcuts. As a is increased further, the re- 
duction of the path lengths by more shortcuts becomes 
less influential, and more shortcuts appear to introduce 
stronger spatial randomness, which eventually makes the 
system less ordered. 

In our simulations, we also notice that the behavior 
of m upon the change of [i is different for small a and 
large a: When a is smaller than some value (a ^ a*), 
the order parameter in the disordered phase m{^ < fj.c) 
approaches zero as N is increased [see Fig. HJa) for the 
WSi network at a = 0.02]. In comparison, at a > a*, 
m{^ ^ fXc) remains finite as N becomes larger as one can 
see in Fig. [D^b) for a = 0.3. It appears that a* is close 
to the value of a at the end point of the lob structure in 
the phase diagram, i.e., a* w 0.1 for the WSi network. 
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FIG. 4: (Color online) Phase transitions in 3D [for (a)-(c)], 4D [for (d)-(f)], and 5D [for (g)-(i)] regular hypercubic lattices in 
terms of the mutation parameter fi. m versus /x, L'^'^m versus ^, and L^''^m versus (^ — fic)L^ are shown in (a), (b), (c) 
for 3D, (d), (e), (f) for 4D, and (g), (h), (i) for 5D, respectively. After the determination of the critical point /Xc in (b), (e), 
(h), all data points collapse to single smooth curve by using the finite-size scaling form in Eq. Q. In 3D, /x^ ~ 5.5, /3 « 0.33, 
and u « 0.63 are obtained in accord with the 3D Ising universality class. In 4D, on the other hand, we obtain fic ~ 5.25, 
/3 « 1/ ~ 0.5, which are the Ising mean-field exponents. In 5D, our simulation results are again consistent with the mean-field 
exponents j3 ^ u ^ 0.5 with fic ~ 6.7. 



At a = 0, WS networks do not possess any shortcuts 
and thus correspond to locally connected regular lattices. 
Observed phase transitions here at a = for both WSi 
and WS2 networks are consistent with the simple expec- 
tation that due to the underlying Z2 symmetry of the 
two defensive alliances, the DAP should belong to the 
same universality class as the equilibrium Ising models 
in ID and 2D regular lattices, i.e., no phase transition at 
finite ^ for ID and the phase transition with the 2D Ising 
critical exponents at finite fj,c for 2D 0, Q- However, the 
existence of discontinuous phase transition at nonzero a 
in WSi and WS2 networks clearly contradicts the above 
simple naive expectation: It has been known that stan- 
dard equilibrium models in statistical mechanics such as 
the Ising and the XY models in the WS networks exhibit 
the mean-field type continuous phase transition, which 
has been attributed to the effective increase of the spatial 
dimensionality due to the shortcuts (see, e.g., Refs. Q)- 
This clearly gives a caveat that one needs to be careful 



in generalizing conclusions drawn for equilibrium models 
to noncquilibrium models. 



The question we are now addressing is if the structural 
inhomogeneity is responsible for the discontinuous transi- 
tion. Since the mean-field theory developed in Ref . |8| can- 
not predict a nontrivial critical point, the nature of the 
transition on the WS network cannot be understood from 
this theory. So it is inevitable to study regular higher- 
dimensional systems, especially beyond the upper criti- 
cal dimensions of the Ising class. For completeness, we 
will study the DAP in regular 3-, 4-, and 5-dimensional 
hypercubic lattices in the next section, before studying 
6-dimensional system in Sec. Ill CI 



B. DAP in 3-, 4-, and S-dimensional regular 
lattices 



The results in three-, four-, and five-dimensional reg- 
ular hypercubic lattices are given in Fig. 3] for hyper- 
cubic lattices with the linear size L (the total number 
of individual species A^ = L"^). We use the system size 
L = 16,20,24,28,32, and 48 for 3D, L = 12,14,16,18, 
and 20 for 4D, and L = 8, 10, 12 and 14 for 5D, respec- 
tively. During numerical simulations, we measure the 
order parameter m in Eq. Q and use the standard finite- 
size scaling form of the magnetization: 



L-^"''f ((m - Mc)^^/'^) 



(3) 



where f{x) is a suitable scaling function with the scal- 
ing variable x, and [3 and v are critical exponents for 
the order parameter and the correlation length, respec- 
tively [lOj . Figure 0] summarizes the numerical results 
for the phase transition in the 3D [for (a)-(c)], 4D [for 
(d)-(f)], and 5D [for (g)-(i)] regular hypercubic lattices. 
Clearly exhibited is the vanishing of the order parameter 
at low ^ [see Fig. |3](a), (d), and (g), in which m versus 
/J, is shown for 3D, 4D, and 5D]. The critical point /Xc is 
determined from the unique crossing point as shown in 
Fig- [11(b), (e). and (h) with L^/'^m versus ^ plotted, and 
then we present the collapse of the numerical data into 
smooth curves in Fig. |4](c), (f), and (i) for 3D, 4D, and 
5D, respectively. In the above finite-size scaling analysis, 
the critical point ^c and the critical exponents are deter- 
mined: ^c ~ 5.5, (3 « 0.33, V « 0.63 in 3D, /.ic ~ 5.25, 
13 « 0.5, V « 0.5 in 4D, ^^ ~ 6.7, [3 « 0.5, v « 0.5 in 5D. 
Accordingly, we conclude that the DAP belongs to the 
same universality class as for the equilibrium Ising model 
in d-dimensions for d = 1, 2, 3, 4, and 5 [l^ (sec Tablc[l|. 



C. DAP in 6-diniensional regular lattice 

We next examine the nature of the phase transition 
in 6D. Due to the practical limitation of the computa- 



TABLE I: Universality classes of the DAP for d-dimensional 
regular lattices, for globally-coupled structure, and for the 
WS network structure. The critical exponents (i and u are 
included when the phase transition is of a continuous nature. 



structure 


tJ.^ 


P 


v 


universality class 


ID regular 


oo 


- 


- 


ID Ising 


2D regular [7, 8] 


6.5 


1/8 


1 


2D Ising 


3D regular 


5.5 


0.33 


0.63 


3D Ising 


4D regular 


5.25 


1/2 


1/2 


equilibrium mean-field 


5D regular 


6.7 


1/2 


1/2 


equilibrium mean-field 


6D regular 


7.4 






discontinuous 


Global coupling [8] 


oo 


- 


- 


no phase transition 


WSi and WSa 


* 


- 


- 


discontinuous 




* ^c for WSi and WS2 depends on the rewiring probability 
a (see Fig. [Sjl. When a = 0, WSi and WS2 are identical 
to ID and 2D regular lattices, respectively. 



FIG. 5: (Color online) (a) m versus /i in 6D. From the size 
dependence of m, we roughly locate /Xc = 7.4(2). (b) Nor- 
malized histogram H{rn) of the magnetization m around the 
critical point for the system size L = 10. The peak position 
suddenly changes around ^c{L = 10) ~ 7.5, indicating the 
discontinuous nature of the phase transition (compare with 
the histogram in Ref. @ for the discontinuous transition in the 
WS network). At large values of fi, the system often exhibits 
asymmetric H(m). Only for convenience of presentation, we 
symmetrized H{rn) to make it an even function of m. For 
comparison, we also show in (c) H(m) for the 5D regular lat- 
tice of the size L — 8. Different from (b), the peak positions 
in H{m) changes smoothly, indicating the continuous nature 
of the phase transition, in accord with the mean-field univer- 
sality class revealed in Fig. 3] 



tional resources, we are limited to use the linear size 
L ~ 7,8,9,10, and 11. Surprisingly, we find that the 
transition nature in 6D is very different from the sim- 
ple expectation of the equilibrium mean-field type and 
becomes discontinuous, similarly to the WSi and WS2 
networks presented in Sec. Ill Al In Fig. OJa), the order 
parameter m is shown as a function of fj., exhibiting the 
transition around fic ~ 7.4. Similarly to the WSi and 
WS2 networks, the change of m near fj,c becomes more 
abrupt as L is increased, which indicates the discontin- 
uous nature of the phase transition. In Fig. El^b), we 
display the normalized histogram H{m) of the order pa- 
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FIG. 6; (Color online) Phase diagram for the DAP in d- 
dimensional regular lattices with the food web given in Fig. [T] 
in the plane of the mutation parameter 1/ ^ and d. Different 
symbols denote different universality classes, and the lines are 
only guides to eyes. For d = 1, 2, 3, 4, and 5, the DAP belongs 
to the Ising universality class. As d becomes larger, the na- 
ture of phase transition is changed to a discontinuous one. 
The point at d = co is included from the result in Ref. [g for 
the globally-coupled case. 



the first upper eritieal dimension dup ~ 4 splitting the 
non-mean-field transition (d < dup ) and the mean- field 

transition (d > dup ). Different from the equilibrium Ising 
model, the DAP model does not always show the mean- 
filed nature for dimensions higher than dup and there 
exists (iii) the second upper critical dimension dup = 6 
beyond which the DAP shows the discontinuous phase 
transition. Unfortunately, we do not have any rigorous 
reasoning for or against the above "conjecture." 

Although the nature of the phase transition cannot be 
determined by the mean-field theory in Ref. [g, it predicts 
that the critical point l//ic in d dimensions should ap- 
proach to from above as d — > cxd, because the globally- 
coupled network is equivalent to the infinite-dimensional 
systems. Consequently, the vanishing value l//ic = in 
both limits of d — > oo and d = 1 makes us expect the 
reentrance behavior as d is increased, which turns out to 
be true as shown in Fig. [51 This is also very different from 
the equilibrium Ising model in which the critical temper- 
ature monotonically increases with the dimensionality. 



III. CLUSTER MEAN-FIELD CALCULATION 



rameter m [m= (cq + C2 + C4) — (ci + C2 + C3) has been 
used to plot H(m)], which clearly exhibits the signature 
of the discontinuous phase transition, similarly to the WS 
network [8|] . As another evidence of the non- mean-field 
nature of the transition in 6D, we also plot L^^'^ as a 
function of fi in the same way as we did for 3D, 4D, and 
5D in Fig.[Hb), (e), and (h)(not shown here). With the 
mean- filed values /3 = i^ = 1/2, the curves do not make a 
unique crossing, which again supports the non-mean-field 
nature of the phase transition in 6D. 

In Fig. [6l we summarize results for the phase transi- 
tions of the DAP in d-dimensional hypercubic regular lat- 
tices. For d < 5, the system shares the critical behavior 
with the Ising model (see Table IJ): No phase transition 
in ID, and continuous phase transitions in 2D, 3D, 4D, 
and 5D with critical exponents corresponding to the Ising 
models in the same dimensions. Very interestingly, as the 
dimensionality d becomes larger (d = 6), the nature of 
the phase transition changes to a discontinuous one. In 
a sharp contrast, the equilibrium Ising model in higher 
dimensions than five belongs to the same mean-field uni- 
versality class as in four dimension (up to logarithmic 
corrections), which defines the upper critical dimension 
of the Ising model {dc — 4). On the other hand, the DAP 
displays very different behavior: Although the equilib- 
rium mean-field universality is identified in 4D and 5D, it 
does not lead to the conclusion that d > 5 should exhibit 
the mean-field universality. From our extensive simula- 
tions of the DAP in the d-dimensional hypercubic regular 
lattices, we propose that the system has three different 
critical dimensions: (i) The usual lower critical dimension 
diow = 1 below which the system is always disordered, (ii) 



From our numerical observation that the discontinu- 
ous transition is not only due to the network property 
but it can also arise from the increased dimensionality 
for hypercubic regular lattices, one would expect that 
there will be a mean-field theory to predict the dis- 
continuous transition. Along this direction, we apply 
the cluster mean-field (CMF) theory [ll| to the two- 
dimensional model. In this section, the lattice point in 
two-dimensional square lattice is denoted by x with the 
decomposition x = niei -1-71262, where n^'s are integers 
and Bi are unit vector along the direction i. The species 
at site X will be denoted by s^. 

The approximation scheme will be detailed in the Ap- 
pendix and here we just sketch the calculation method. 
To begin with, we write down the exact evolution 
equation for the marginal probability (see Appendix) 
Piisx, Sx+ei, Sx+e2, Sx+ei+ei) ^om thc mastcr equation. 
Since each s can take six different values (s = 0, 1, • • • ,5) 
we have to find out 6^ equations. Thanks to the spa- 
tial (rotation and inversion) and intra-alliance symme- 
try [-P({s}) = P{{s + 2}) with modulo 6], the number 
of equations is greatly reduced to 77, and the probabil- 
ity conservation condition reduces it to 76. To get the 
steady state magnetization, we numerically iterate those 
76 equations until the magnetization does not change sig- 
nificantly. The results are sunrmarized in Fig. [T] We then 
use the fitting function ?ti(/.j) ~ (l//^c — Vm) with the 
mean-field critical exponent (3 = 1/2 and estimate the 
critical point Hc — 6.41, which should be compared to 
the Monte Carlo simulation results fic — 6.50(4) in Ta- 
ble [l] and Ref. [^. Unlike the higher-dimensional systems 
with d > 6, the CMF does not seem to predict the dis- 
continuous transition although it gives us an accurate 
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FIG. 7: (Color online) Magnetization m as a function of 
1/fi in the cluster mean-field approximation. In the fitting, 
m(/Lt) = A(l/fj.c — 1/ ^Y'"^ is used with two fitting parameter 
A and ^c, resulting in /ic — 6.41. 



estimation of the critical point for the 2D system. How- 
ever, CMF does not provide a clue to why the models 
in higher dimensions above 5 as well as in WS networks 
display the discontinuous transition. It is quite difficult 
to believe that the CMF with a larger cluster size than 
2 can provide an explanation, because the actual transi- 
tion in two dimensions is still continuous. One way might 
be applying the CMF to 6-dimensional systems, but it is 
simply too complicated to calculate. 



has been observed as the rewiring probability a is in- 
creased at a fixed mutation rate, implying the intricate 
role of the random shortcuts on the phase transition. 

Hypercubic regular lattice structures in d dimensions 
have also been used as underlying spatial interaction 
structure of the DAP. Observed is that the phase dia- 
gram in the plane of the dimensionality and the mutation 
rate again displays an interesting reentrant behavior (sec 
Fig- ©7 which is in a striking contrast to the equilibrium 
Ising model. For the latter model, the critical tempera- 
ture is simply an increasing function of the dimensional- 
ity. We have identified the universality class for various 
dimensions {d = 1, 2, • • • ,6) and d = 1, 2, • ■ • ,5 exhibits 
the same critical behavior as for the equilibrium Ising 
model. In contrast, as d is increased further, discontinu- 
ous phase transition has been observed for d = 6. In the 
hope to find a theory predicting the discontinuous transi- 
tion, we have applied the cluster mean-field approxima- 
tion for the two-dimensional systems, but we only find 
the continuous transition within this scheme. Still, the 
reason why the higher-dimensional systems exhibit the 
discontinuous transition remains mystery, which can be 
an interesting theoretical question to be pursued further. 
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IV. SUMMARY 

We have investigated the alliance breaking phase tran- 
sitions of the six-species predator prey model called the 
defensive alliance process (DAP) in the complex networks 
and also in the d-dimcnsional regular lattices. Identified 
nature of the phase transition is summarized in Table U 
with critical exponents included when available. 

In the WS network, we have observed that the alliance 
breaking phase transition is of a discontinuous nature. 
Interesting reentrant phase transition in the WS network 



APPENDIX: TWO-DIMENSIONAL CLUSTER 
MEAN-FIELD THEORY 

In this Appendix, the equations and symmetry of the 
cluster mean-field approximation with the cluster size 
two will be detailed. We start from deriving the exact 
equation for the marginal probability for the local con- 
figuration at four sites x, x + ei, a; + 62, and a; + ei + 62 
(x is the vector designating a lattice point and e^'s are 
unit vectors along i-th direction). 

The marginal probability P4 is defined by 



J 



P4(a, 6, c, d) = P4 



a b 
c d 



-E' 



p( 



3i ^j ^x-\-e2 ^7 '^x-l 



= 6, 



(A.l) 



r 



where the time dependence is implicitly assumed and 
the primcd-sum (^ ) run over the all possible configura- 
tions with four specified sites having the denoted species 
(0 < a,b,c,d < 5). We always assume the modulo-6 
equivalence among species indices, e.g., species 6 is equal 



to species 0, and so on. Due to the translational invari- 
ance, P4 does not depend on x if initial condition does 
not have x dependence or if we are only interested in the 
stationary state. Hence we can safely omit the explicit x 
dependence for the function P4. 



The exact time evolution for P4 can be written as 
2 
dt'" \c d) ~ ^ 

t — 

a 




(Sab + 5ac)PA ( ^ ^ J + {5ab + 5m)Pa ( ^ J + {^ac + 5cd)PA L + i d 



+(^M + Mn(;,^^) 



y^(l - 5iz){5a.b+i + 5a.c+i + Sb,d+i + Sc,d+i)P4 i ^ ^ 



a b^ 



(A.2) 
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In Eq. (jA.2[) . we also introduced P5 which is the marginal 
probability of the local configuration with 5 sites, taking 
the form given as an argument. Due to the rotational 
and mirror symmetry, one can easily see that 




d\ ^P, 



e c a 
d b 



= A 



d b 
c a 



(A.3) 



Hence we do not have to deal with 8 different local con- 
figurations independently while treating P5. 

Due to the hierarchy appearing in Eq. (|A.2[) (P5 is not 
reducible in terms of P4), it is not easy to solve it exactly, 
which necessitates the use of the approximation scheme. 
To this end, we treat the correlation within squares of 
linear size 2 completely and neglect the correlation be- 
yond linear length 2. To be specific, we are using the 
approximation scheme such that |ll| 




Pi 
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c d 
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Pi 



d) 



(A.4) 



where P2 and P3 are marginal probabilities defined simi- 
larly to P4. The probability conservation makes it possi- 
ble to find P2 and P3 once we know P4 from the relations 



1=0 



e I 



p,{cd) = j2p. 



c d 



(A.5) 



(A.6) 



The above approximation scheme along with the symme- 
try consideration given in Eq. (jA.3|) lets Eq. (jA.2p closed 



in the sense that the number of equations are equal to 
that of variables. Since it is still infeasible to solve the 
approximate equation analytically, we resort to the nu- 
merical solutions. 

One may treat the 6^ = 1296 equations directly with- 
out considering the degeneracy of Pj's. However, sym- 
metry consideration reduces the number of equations we 
have to deal with considerably. Now we will show that 
actually we have only to treat 76 equations to get the full 
solution. 

There are three symmetry operations which are sum- 
marized as follows: 
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c d 

a b 
c d 

Pi 



Pi 



Pi 



b d 
a c 
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Pi 
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d b 
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b d 
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c d 



= Pi 



a + 2 b- 
c+2 d 



(A.7) 
(A.8) 
(A.9) 



which are rotation (±90°), mirror (against the diagonal 
axes), and intra-ally symmetric operations, respectively. 
To find the degeneracy due to symmetry, let us first cat- 
egorize the local configurations according to the number 
of same species among 4 sites. There are 5 such cat- 
egories which take the form {0,0,0,0} (6), {0,0,0,1} 
(120), {0,0,1,1} (90), {0,0,1,2} (720), and {0,1,2,3} 
(360); where the species numbers in curly brackets are 
just for the representative purpose, the order of four ele- 
ments in the curly braces are irrelevant, and the numbers 
in parentheses indicate the total number of local configu- 
rations of the corresponding categories whose sum should 
be64. 

By considering the above symmetry operations, we 
reduce the number of independent variables. When 



a = b = c = d, the intra-ally symmetry reduces the inde- 
pendent variables from 6 to 2. When three of four species 
are the same such as a = b = c ^ d, rotational sym- 
metry as well as the intra-ally one reduces the number 
of independent variables. For example, P4{a,a,a,d) = 
P4^{a,a,d,a) = Pii{a,d,a,a) = Pi(d,a,a,a) by rotation 
and P4(a, a^a^d) = ^4(0 + 2, a -I- 2, a -I- 2, d -I- 2) by intra- 
ally transformations. From this consideration, one can 
find that there are 10 independent variables out of 120 
variables. In case a = b ^ c= d, the symmetry consider- 
ation shows that (we omit P4 for convenience) 



a a 
c c 



a c 
a c 



c c 

a a 

a c 

c a 



c a 

c a 

c a 

a c 



(A.IO) 



which reduces the independent variables from 90 to 30. 
The intra-ally symmetry once again reduces the number 
from 30 to 10. 

The fourth case is a = b ^ c ^ d ^ a. The symmetry 
enforces the equivalence among configurations such that 



a a 
c d 



a d 
a c 



d c 

a a 

a c 

d a 



c a 

d a 
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(A.ll) 



where exchanging c with d, which is equivalent to the 
mirror transformation, also gives the equivalent configu- 
rations. Hence Eq. (|A.11[) combined with the intra-ally 
symmetry reduces the number from 720 to 40. The last 
category contains all different species. This category has 
(^) = 15 elements. Since the rotation and the mirror 
symmetry operations always conserve the diagonal rela- 
tions, each set has 3 different classes. Due to the intra- 
ally symmetry, however, only 15 configurations remain 
independent. To sum all independent configurations, we 
find that there are 77 (=2-1-10-1-10-1-40-1-15) independent 
variables. Due to the probability conservation, one equa- 
tions becomes reducible by other variables. So the final 
number of independent variables is 76. 



One may solve the stationary state solution by setting 
d/dt = in Eq. (jA.2[) . For us, this turned out not to 
be easy, so we numerically integrate Eq. (|A.2p starting 
from the fully magnetized state with intra-ally symmetry 
such that P4(a, 6, c, d) = 1/3" with n to be the number 
of members in alliance I among a, b, c, d and found the 
stationary state magnetization, which is summarized in 
Fig. [71 
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